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FOREWORD 


This  support  analysis  and  scientific  programming  was  per- 
formed under  contract  to  the  Analysis  and  Simulation  Branch 
( SUA ) of  the  Air  Force  Geophysics  Laboratory  (AFGL), 
Hanscom  Air  Force  Base,  Massachusetts.  Programs  described 
herein  and  documentation  of  such  can  be  obtained  from  the  com- 
puter library  of  SUA. 
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SHARP  CUTOFF  NON-RECURSIVE  FILTER  DESIGN 


Initiator 
Problem  No. 
Project  No. 


A.  Slobodnik,  Jr. 

3057 

5635 


The  purpose  of  this  project  is  to  design  a sharp  cutoff,  flat  top,  bandpass 
filters  using  a Lagrange  interpolation  design  algorithm.  Several  articles  on 
this  subject  were  gathered  and  read.  Possible  computationally  efficient  in- 
terpolation schemes  were  examined. 

The  technique  to  be  investigated  initially  uses  trigonometric  polynomial 
approximation.  Given  desired  values  for 

N = order  of  polynomial 
Np  = maximum  number  of  passband  ripples 
Sp  = maximum  deviation  allowable  in  the  passbr 
Ss  = maximum  deviation  allowable  in  the  stopband. 

The  design  algorithm  first  estimates  an  initial  set  of  extremum  frequen- 
cies F^\  K = 1, ...  ,N  + 1,  at  which  the  extrema  of  the  desired  frequency 
response  Hn(F)  are  to  be  located.  It  then  uses  the  barycenter  form  of  the 

v x-L. 

Lagrange  Interpolation  formula  to  obtain  an  N—  - order  polynomial  that 
passes  through  the  maximum  ripple  values  (1  +Sp  in  the  passband  and  +Ss 
in  the  stopband)  at  these  prescribed  frequencies.  In  general,  this  first 
Lagrange  polynomial  will  achieve  its  extrema  at  frequencies  distinct  from 
the  estimated  set  of  original  frequencies  F^  and  additional  stages  of  the 
algorithm  will  be  required. 

The  final  algorithm  is  a slight  extension  of  an  algorithm  due  to  Parks  and 
McClellan.  Basically,  a weighted  Chebyshev  filter  with  trigonometric  co- 
sine polynomials  is  found  for  the  step  function, 


11(f)  *•-  1 , 0 £ f £ FP 

H(f)  = 0 , FP  + DELF  £f  £l/2  , 

where  DELF  < 0 is  introduced  by  the  algorithm  and  adjusted  iteratively  until 
the  R8,  RB  tolerances  are  satisfied.  Excellent  results  have  been  obtained 
with  this  approach. 


The  results  were  submitted  to  the  investigators 


ANALYSIS  OF  ELECTRIC  FIELDS  IN  VISUAL  AURORA 


Initiator 
Problem  No. 
Project  No. 


M.  Smiddy. 

4535 

8617 


This  problem  deals  with  the  determination  of  the  ambient  electric  field 
vector  as  measured  during  a series  of  rocket  flights  into  a visual  aurora. 

In  the  geomagnetic  polar  region,  where  the  earth's  magnetic  field  (carry- 
ing trapped  radiation  belt  particles)  dips  into  the  atmosphere,  plasma  theory 
dictates  that  an  electric  field  should  be  created.  The  electric  fields  should 
be  very  weak  - a few  tens  of  millivolts  per  meter. 

The  low  values  of  the  electric  fields  make  measurements  difficult  and  the 
exact  locations  of  the  field  themselves  are  uncertain.  Visible  auroras,  how- 
ever, give  clues  as  to  where  to  lode  for  them. 

A Black  Brant  VC  rocket  was  launched  on  April  4,  1972  into  a visible 
aurora  display  over  Fort  Churchill,  Canada. 

Detection  of  the  fields  were  made  by  six  solid  spheres  mounted  on  three 
pairs  of  diametrically  opposed  booms.  The  spheres  reach  their  respective 
local  plasma  potential,  due  to  the  high  plasma  conductivity.  The  potential 
difference  between  any  sensor  pair  gives  the  electric  field  component  along 
the  boom  direction.  A spin  was  imparted  to  the  vehicle  about  its  longitudinal 
axis  to  enable  the  electric  field  to  be  measured  anywhere  in  a direction  nor- 
mal to  this  axis.  Due  to  the  slow  precession  of  the  longitudinal  vehicle  axis 
itself  the  plane  in  which  the  electric  field  is  sampled  changes  in  space. 

Each  pair  of  booms  will  be  analyzed  separately.  The  first  stage  of  the 
analysis  consist  in  processing  and  merging  the  magnetic  tapes  containing  the 
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attitude  data  of  the  rocket  and  the  booms,  and  the  tapes  containing  the  values 
for  the  potential  differences  for  each  sensor  pair.  These  merged  tapes  re- 
present the  input  to  the  program  written  to  analyze  the  ambient  magnetic 
field. 

Following  is  a functional  description  of  the  computer  program  written  to 
analyze  the  data  and  a description  of  the  mathematical  method  used. 

The  ambient  electric  field  2A  exists  in  space  as  a function  of  position  and 
time,  E.  = lL(r,t).  Probes  mounted  on  booms  extending  from  a rocket 
measure  a potential  difference  AV,  when  flying  through  regions  of  the  iono- 
sphere in  which  an  electric  field  E exists.  If  the  probes  are  aligned  along 
an  axis  H(t).  moving  with  velocity  v(t),  through  the  earth’s  magnetic  field 
B(r),  then  the  total  voltage  measured  is 

AV*S*S=Sa*S+vx3*S  . 

A 

The  object  of  the  program  is  to  reconstruct  1!A  from  the  experimental 
data  consisting  of  AV,  S/|R|,  v,  and  B,  all  versus  time,  v and  § are 
obtained  analytically. 

The  components  Ex,  Ey,  and  Ez  of  the  ambient  electric  field  can  be 
approximated  from 

2a.^=AV-vxS*S 

A 

and  the  direction  cosine  of  &(t)  by  fitting  in  the  least-square  sense  to  a 
polynomial  of  arbitrary  degree  n.  The  fitted  field  ®A  can  be  output  as  a 
list  of  the  components  vs.  time  or  in  form  of  Calcomp  plots  of  the  compo- 
nents vs.  time. 

The  program  is  extensively  documented  and  is  logically  broken  down  into 
the  following  parts: 
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I.  Initialization  phase.  Reads  various  parameters  and  initializes  for  each 
run. 

II.  Read  R data.  Reads  attitude  data  from  tape  and  sets  up  the  array  of 
time  points. 

III.  Read  E data. 

IV.  Least-square  analysis.  Perform  the  least-square  analysis. 

V.  Output  phase.  Sets  up  printed  output  and  writes  required  plots  on 
plotter  tape. 

The  electric  field  2(t)  is  defined  as 
3 

E(t)  Ej(t)  jTj  (CMj}  o.n.  basis  of  R)  . 

i=l 

The  object  is  to  determine  the  components  E^t)  given 
y = 2 • ft  and  3(t)  . 


We  approximate  E^(t)  by  a polynomial  in  t 

Pi 

V'-IX1  • 

Define 


n r 3 


-i2 


-E  2>y  E‘‘ 

Pi  i=i 


■i  ‘ y(V 


by  (l) 


N 


S = 


k=l 


3 ^ A 

r 

k=0  / 


-.2 


y(tj) 


(i) 


where  N = Number  of  time  points.  Then  calculate 


SS  _n 


(least  squares) 


which  will  produce  M = Ej*=1(Pj  + 1)>  equutlonB 


The  coefficient  of  is 


L’Wj'vy'?  • 

j=l 

These  coefficients  are  the  elements  of  an  M x M matrix  which  together  with 
the  M right-hand  sides  form  a system  of  M linear  equations.  This  system  is 
solved  for  file  M unknowns  by  the  method  of  matrix  inversion.  The  SSP 
routine  MINV  is  used  to  achieve  this. 


The  rocket  Involved  was  spinning  around  Its  axis  at  a much  faster  rate  than 
expected  and  this  made  it  practically  impossible  to  extract  useful  data  for  the 
Intended  analysis.  It  was,  therefore,  decided  by  die  investigator  to  terminate 
the  effort  and  analyze  another  rocket  flight,  reported  under  Problem  No.  4764. 


NET  BP  AC  K-CHARGE- LIMITED  CURRENT  FLOW 


Initiator  : P.  Gianino. 

Problem  No.  : 4630 

Project  No.  : 6621 


This  problem  is  aimed  at  solving  nonlinear  Vol terra  integral  equations 
of  the  first  and  second  kind,  numerically,  for  the  net  space-charge-limited 
current  flowing  between  2 plane-parallel  plates  emitting  secondary  electrons. 
These  electrons  are  emitted  when  radiation  impinges  in  the  metal  plates. 

The  problem  consists  of  5 parts:  (a)  compute  an  upper  threshold,  (b) 
compute  a lower  threshold,  (c)  compute  a threshold  for  20  - 40  points  be- 
tween the  upper  and  lower  threshold  (the  space-charge-limited  region), 

(d)  (e)  evaluate  some  expressions  which  depend  on  the  answers  obtained 

in  (a)  and  (b). 

A typical  integral  that  occurs  in  (a),  (b),  and  (c)  is  of  the  form 
n2  . J_  f'402VE  d» 

415  J»  ^*,u)  + Le-402VEh»  ’ 

where  hA(u)  = eu  - 1 + 2«/u/ff  + eu  erf«/u,  L and  D are  parameters  with  known 
values.  The  problem  is  to  find  the  value  of  VE  for  which  equality  holds,  i.  e. , 
VE  is  the  desired  threshold. 

There  are  several  problems  which  must  be  overcome  before  the  problem 
can  be  solved: 

(1)  the  integrand  has  a singularity  at  u = 0; 

(2)  the  problem  is  highly  nonlinear; 


(3)  each  change  in  VE  changes  the  integrand; 

(4)  the  expression  eu  erf</u  in  h^u)  creates  substantial  numerical 
problems. 

(2)  and  (3)  are  surmountable  since  the  integrand  and  the  Integral  are  both 
monotonic  functions  of  VE.  This  allows  a use  of  some  well-known  iterative 
techniques. 

(1)  was  overcome  after  a substantial  study  of  the  sensitivity  of  values  of 
the  integral  to  changes  in  the  lower  limit  (i.e, , replacing  0 by  c and  finding 
an  c suitable  for  the  purposes  of  this  study). 

(4)  was  the  most  difficult  to  get  around.  A number  of  published  formulae 

o 

were  used  with  little  success.  A continued  fraction  expansion  for  e*  erf  X 
was  coded  and  discovered  to  be  wrong  for  X « 1.  Finally,  a combination  of 
published  and  unpublished  results  was  coded  and  incorporated  successfully 
into  the  program. 

In  solving  parts  (1)  and  (2),  the  problem  was  treated  in  two  parts: 

(a)  Find  2 values  of  VE  for  which  the  right-hand  side  was  greater 

o 

and  less  than  the  desired  value,  O . 

(b)  Apply  the  technique  of  "interval  halving"  to  the  Interval  found 
in  (a),  to  obtain  the  solution  VE. 

A similar  approach  was  used  in  (c).  All  integrations  were  performed  by 
using  a 32  point  Gauss-quadrature. 

The  problem  was  completed  and  the  results  given  to  the  investigator. 
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ELECTRON  CAPTURE  FROM  TI+  AND  Li* 


Initiator 
Problem  No. 
Project  No. 


R.  Mapleton 

4666 

8627 


THE  PHYSICAL  PROBLEM 

This  problem  is  concerned  with  estimating  the  cross  section  tor  the  reaction 

H(n)+^Tl++(m)  (1) 

n m 

for  energies  in  the  range  1 < E s 16  keV,  where  n denotes  the  states  of  atomic 
hydrogen  and  m denotes  the  states  of  the  doubly-ionized  Thallium  atom  Tl++ . 
Since  capture  into  H(1S)  dominates,  we  can  write  H(n)  = H(1S)  in  (1);  also, 
since  for  these  impact  energies,  cross  section  for  capture  of  the  outermost 
electron  dominates,  we  can  simplify  (1)  to  obtain 

H*  + Tl+  (6S2  : XS)  — > H(1S)  + Tl++  (6S : 2S)  (2) 

using,  respectively,  1 electron  and  2 electron  ions  for  Tl+  and  Tl++. 

We  next  construct  a Schrodinger  wave  equation  for  (2).  For  the  impact 
version  of  this  problem,  the  equation,  in  atomic  units,  is 

(h-1  £)*>  = <>  , (3) 

where  H denotes  the  Hamiltonians.  After  capture,  the  Hamiltonian  is  given  by 

j j V^*l) 

H ■ H(pl)  + H(C**2)  - (4) 

r2p  r12 


H*  + T1+(6S2  :XS) 
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and  <"<t<iiiy  (before  capture)  the  Hamiltonian  H la  given  by 


H = H(C+ 12)  -~r — , (5) 

r2p  'lp 

where  H is  the  Hamiltonian  for  the  hydrogen  atom,  the  T1++(6S),  the  doubly 
ionized  Thallium  atom  and  T1V+(6S2),  the  singly-ionized  Thallium  atom. 
Also,  in  (5),  p denotes  proton,  subscript  1 denotes  electron  1,  subscript  2 
denotes  electron  2,  C++  denotes  the  Tl++  core  and  C+  the  H+  core.  Dis- 
tances between  particles  are  denoted  by  r12,  rpl,  rp2,  rln,  and  r^  with 
n denoting  the  T1  nucleus. 


The  Hamiltonians  In  (4)  are  given  by 


and 


where  lVf(C++2)l/2  is  derived  from  the  eigen  equation 
H(C4"*,2)  $(T1++)  = -€(++)  ftTl++)  . 


(6) 


(7) 


(8) 


We  now  change  coordinates,  from  the  hydrogen  atom  to  the  Tl+  electron 
coordinates  r^,  and  obtain  from  (4): 


H 


( M IJ  V^l)  V,(C~2)  \ l 

■ i ■$—  2 +rj  rlp  r ^ 


(9) 


i 

£l2"Enl-£n2  \ 

/ 


J = 1,2  5 rnJ  ~ rJn  * (10) 
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Apart  from  the  Coulomb  Interactions,  (9)  represents  the  Hamiltonian  for 
the  function  *2(Tl4  4«S>.  Replacing  Vj(C++j)  by  Vt(C+J)  - ( l/r lg)  in  (9),  we 
have 


a.(.i+  .1+ 

\ 2 " 2 112  2 2 r12/  rlp  r2p 


(U) 


"Vs 


rlp  r2p 
and  we  approximate  Hj  by 

Hi^(Tl+)  = -e(+)  ^Tl+)  , 


(12) 


(13) 


where  *(+)  denotes  the  (experimental)  energy  of  T1  . Next,  the  assumption 
is  made  that  the  interaction  of  e(H)  with  the  Tl++  ion  is  represented  by 
Vj(C+l),  imagining  T1+(6S2)  as  constructed  of  T1++(6S)  and  e(H),  we  get 

Vf(C++j)  V (C+j)  . 

—2 = -4 f-  , J=M  , 

2 2 r12 

which  can  be  seen  to  be  exact  at  large  distances  from  the  nucleus. 


(14) 


We  now  consider  the  scattering  problem:  Let 


*=  A(t)  ^ + B(t)  tf>f 


(15) 


where  ^ and  ^ represent  the  wave  functions  of  the  system  before  and  after 
scattering,  respectively.  These  are  given  by 

*t  * X*JlB'£2n>  “tp{1  [(-*  ' t)  * - T (5l  + 52>|  <“> 

*t  • x(Ipi'  W “*•  {*  [(f  1 + ' M ' T") 1 * i <5l  ' 52,]|  - <l7> 


and 


13 


where 


The  scattering  equations,  from  (3)  are  then  given  by 


From  (11),  we  have 


and  similarly 


Lis  .is  =0 

1 2yls  2^8  r12  2 2 l*)vi  Qa 


(24) 


Now,  writing 


and 


(25) 

(26) 


we  then  can  obtain,  from  (15)  and  using  (18)  - (24),  the  following  coupled 
system  of  (complex)  ordinary  differential  equations: 

1(1  - |S|2)  A(t)  = [Vu  - S*Vfl)  A(t)  + (Ulf  - S*Uff]  B(t) 

i(l  - |S|2)  fi(t)  = [Vfl  - SVjj]  A(t)  + [Uff  - SUtf]  B(t) 

with  the  initial  conditions 


(27) 


where 


and 


A (t  = -»)  = 1 1 
B (t  = -»)  = 0 f , 

V,(C*1)  t l v^c^)  v^c^) 


r12  r2p  2 


and  we  write 


Vu'/070"*’V*IdrInd'': 

v«7„T*f*v**dri"d^ 


(28) 


(29) 


2n 
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(Continued) 


*i  rf  In  2n 


(30) 


CD  « CO 

u“i  I 


and 


S'/0  J[  *?  *1  drla  dr2n  • 

Thus,  to  obtain  the  quantity  of  interest,  the  cross  section  for  capture  from 
T1+(6S2)  into  H(1S)  + Tl^Ss)  by  protons,  we  have  to  evaluate 


/•OO 

Q = 4irp2Jo  p|B  (t  = *)|2dp  , (31) 

where  R2  = p2  + v2t2,  v = 0.2</E  and  t is  time,  p is  the  impact  parameter, 
and  v is  the  initial  velocity  with  respect  to  the  Tl+  ion  at  rest  and  E is  the 
energy.  This,  in  turn,  requires  the  solution  of  the  initial  value  problem, 
given  by  (27),  (28)  for  sufficient  values  of  the  Impact  parameter  p to  enable 
the  integration  in  (31)  to  be  carried  out.  Note  that,  in  turn,  the  solution  of 
(27),  (28)  is  in  itself  extremely  difficult  due  to  the  complexity  of  the  elements 
of  the  coefficient  matrix,  given  by  (30).  Thus,  to  make  the  problem  at  all 
numerically  tractable,  a number  of  simplifications  must  be  imposed. 


+ ++ 

First,  since  the  radial  wave  functions  for  T1  and  T1  and  the  interac- 
tions Vj(C+j)  and  V^C^j),  j * 1,2  are  given  in  numerical  form,  it  is  not 
possible  to  simplify  the  matrix  elements.  However,  if  momentum  transfer 
is  neglected  (i.  e. , set  to  zero)  then  we  can  obtain  the  simplification  given 
by  omitting  the  factors  exp  [(-i/2)  vfS^  + S2)]  and  exp  [(-i/2)  v(Sj  - Sg)] 
from  the  expressions  for  (|)j  and  >Pp  respectively,  in  (16)  and  (17),  yielding 


*f  ' X(Ipl>Ip2>  exp  [‘  («1  +»*  * t)  *]  ’ 


(32) 
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Next,  since  the  integrals  In  (30)  are  effectively  double  Integrals,  we  will 
expand  the  Integrands  in  Legendre  series,  using 

rpj=rnj”R  O3) 

2 2 2 2 

with  R = p + v t , etc.  The  independent  variables  (rln,r2n,t)  are  used 
since  the  wave  functions  are  given  numerically  in  these  variables.  The 
Legendre's  expansions  used  are 


|r  -r  r1=-A_ 
1 nl  n2 1 rni> 

exp(-|rni-R|)  j 


lrnl-Rl 


2r  ,R 


(er<  - e_r<)  e“r> 


er*-e-r</  , l\  /e-*  + e-r<\  -r> 
[ *>)~\ *> )e  » 


where  < and  > denote,  respectively,  the  lesser  and  greater  of  the  two  vari- 
ables on  the  left-hand  side.  The  Auctions  P(T1+),  tp(Tl ++)  correspond  to 
rR(r),  where  R(r)  is  the  corresponding  radial  function.  Also,  the  functions 
P and  >P  are  normalized,  i.  e. , 


f |P2(x)|dx=f  |^2(x)  | dx  = 1 , 

Jo  Jo 


and,  as  well  as  the  interactions  V^,  are  given  numerically.  Then,  using 
(34)  in  (30),  the  matrix  elements  become 


V-e-y 
. y 1 
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mamm. 


vfl  - -S1(R){1  +k/0R  *<x>  p«  *} 

!R  -R 

pw e * 


-R  rR  „ 

v(  <**- 


e“xl  P(x)  dx  } I 4>(y)  P(y)  dy  . (38) 


uif  * { /R  p<x>  x [e  'r  H)'  *“]  e'X 


' 2R  -2R 
e - e 


' 2x  -2x 
e - e 


... 


(Continued) 
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“ Jr  ^y^dy  (40) 

where  x and  y are  dummy  variables  of  integration. 

We  next  separate  A and  B into  the  real  and  complex  parts  with  the  result 
that  we  now  have,  a coupled  system  of  four  ordinary  (real)  differential 
equations  in  four  unknowns.  We  have 


A = C + D.  , B = E + Fj  , (41) 

hence,  (27)  and  (28)  become  the  initial  value  problem 

(1  - |S I2)  6(t)  = I11D(t)  + I12[-E(t)  siny  + F(t)  cos  y] 

-(1  - Is  |2)  £>(t)  = l„C(t)  + I19[E(t)  cos  y + F(t)  siny] 

o (42) 

(1  - |S  f ) E(t)  = I21[C(t)  sin  y + D(t)  cos  y]  + I22F(t) 

-(1  - |S  |2)  f'(t)  = I21[C(t)  cos  y - D(t)  sin  y]  + Ig2E(t) 
with  initial  conditions 

C(— o»)  = 1 ; D(-®)  = E(-®)  = F(-»)  = 0 . (43) 


Note  that  the  variable  y arises  from  (16)  and  (17)  and  is  defined  by  (20),  i.e. , 


(44) 


Also,  that  we  have 


rn  = -vu + SVA  - suit + 
X12  = X21  = Ufi  " SUff 
r22  = Uff 
Ufi  = Uif 


(45) 
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with  V14>  S,  etc.  given  by  (40). 


Finally,  note  that  it  can  be  shown  that 

T(t)  = |C(t)  |2  + |D(t)  |2  + |E(t)  |2  + |F(t)  |2 
+ 2S{[C(t)  E(t)  + D(t)  F(t)]  cos  y 
+ [C(t)  F(t)  - D(t)  E(t)l  siny} 


f 

I 


-• 

f 


evaluated  at  t = •,  is  equal  to 


i.  e. , at  t = 0,  for  the  case  of  momentum  transfer  included,  we  have  T(<*>)  = 1, 
hence,  the  closeness  of  the  calculated  value  of  T at  t = <*>  to  1,  indicates  the 
correctness  of  the  numerical  solution,  with  the  simplifications  as  detailed 
above. 

Programs  were  developed  to  solve  these  equations  and  the  results  were 
presented  to  the  investigator. 


STUDY  OP  A DIELECTRIC  SLAB  IN  A WAVEGUIDE 


Initiator 
Problem  No. 
Project  No. 


H.  Steyskal. 

4669 

4600 


This  problem  studies  the  mode  propagation  for  a waveguide  loaded  with 
a dielectric  slab. 


A program  was  written  which  iteratively  obtains  a normalized  complex 
wave  number,  0 , which  is  determined  from  the  following  equation 

«(n)  «(n-l)  + ^<n-l) 

m m m 

for  each  mode  m,  where 

m = 1, 2, ... , M as  an  outer  index 
n = 1, 2, ... , N as  an  inner  index  , 

and  where  (i^  is  the  complex  conjugate  of 

fe~-\ (if  • 

If  m = even  number,  then 


B = cos  (qc) 


SiSLfifc-lSU  [loos  (pb)  oo.  (pa, -qSSLm  •fiSiESi.j  , 


If  m = even  or  odd,  then 


Note  that  everywhere  the  following  are  complex  numbers:  p,  0,  A fi,  D, 
Dc»  D0>  B»  Bc»  B/J* 

The  results  were  given  to  the  investigator, 

I 

I 


I 


ANALYSIS  OF  A WAVEGUIDE  JUNCTION 


Initiator  R.  Mailloux 

Problem  No.  : 4670 

Project  No.  : 4600 


The  purpose  of  this  problem  la  to  produce  a modal  analysis  of  a waveguide 
Junction.  The  program  solves  the  problem  of  the  junction  of  three  waveguides. 
Each  waveguide  is  partially  loaded  with  dielectric  aa  shown  below. 


The  parameters  M and  N are  related  to  the  number  of  modes  in  die  large 
waveguide  and  the  small  waveguides.  Aptitudes  and  aj  are  the  small  wave- 
guide excitation  coefficients  and  r is  the  relative  dielectric  constant  of  the 
dielectric  slabs. 

A system  of  N linear  equations  (corresponding  to  the  N modes  in  the  wave- 
guide) are  solved. 
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The  above  values  of  the  B equation  are  various  input  pa  rameters  or  functions 
of  input  parameters. 

\ 

Work  is  continuing  on  this  job.  The  program  written  is  now  being  put  through 
various  production  runs  for  Mr.  Robert  Mailloux  who  is  taking  over  for  Mr. 

Hans  Steyskal. 

We  then  combine  the  scattering  matrix  data  produced  from  this  (the  "junc- 
tion" program)  with  some  data  of  the  initiator  to  arrive  at  the  reflected  Helds 
for  the  "junction”  in  the  presence  of  a phased  array. 

A program  was  developed  to  perform  this  and  the  results  were  given  to  the 
investigator. 


ROCKET  EXHAUST  PLUME  STUDY 


Initiator 
Problem  No.  : 
Project  No.  : 


R.  E.  Good  LKC 

4683 

IUR 


The  definition  of  high  altitude  rocket  exhaust  plumes  has  been  developed 
by  numerous  investigators  and  published  in  various  sources.  It  remains 
to  solve  numerically  the  formulation  given  by  Hill  and  Jarvinen.  The  defi- 
nition of  the  plume  shape  will  be  used  in  the  determination  of  missile  thrust. 


The  problem,  then,  is  to  solve  numerically  the  following  equations: 
1(0,9)  = 0 exp  (-A2!!2)  sin  8 cos  (8  - <p)  d3 
0 = 1-  cos  (0  + a)  , A = * ’ (^X1  + 

L - a[>  + < V » t]  H - » ”3“e]'1/2 


F=  [2/1  y=  + 3 to2 

L l2y-ti  ry.) 


L exp  (-A2B2)  sin2 


(0-9) 


sin  0 


R -flL 

Rc  F * 


The  inputs  are  D/T,  Mg,  Y Y^,  M^,  and  a.  The  plume  shape  is  derived 
using  the  definition  of  the  radius  of  curvature  in  Cartesian  coordinates.  The 
local  slope  of  the  contact  surface,  tan  9,  yields  the  first  derivative.  The 
radius  of  curvature,  when  solved  by  difference  method,  yields  the  change  in 
slope: 


By  knowing  p initially,  R,  is  computed.  From  Rq,  ¥ Is  determined.  This 
permits  definition  of  <p  to  start  the  next  integration  step. 

The  program  was  developed  and  the  results  presented  to  the  investigator. 


INTERPRETATIONS  OF  8P ECTROSCOPIC  DATA 


Initiator  : J.  Mansan 

Problem  No.  : 4715 

Project  No.  : 6688 


This  problem  is  one  of  explaining  experimental  data  from  a spectroscope 
in  a rocket.  Photon  counts  should  Increase  at  or  near  known  laboratory 
spectral  lines,  but  can  vary  in  intensity;  be  contaminated  by  other  lines 
nearby;  suffer  from  general  background  noise;  change  due  to  noise  bursts, 
mechanical  difficulties  over  parts  of  a scan,  or  even  variations  of  the 
transfer  Amotion  of  the  aperture  during  the  flight;  or  any  combination  of 
the  above. 

Since  the  problem  involved  human  judgments  about  graphical  data,  it  be- 
came a prime  candidate  for  solution  via  interactive  graphics.  A program 
was  written  which  plots  the  experimental  points  on  a CRT  screen.  The  user 
then  inserts  spectral  lines,  and  is  allowed  to  vary  their  intensity,  shift 
their  wavelength,  and  add  or  subtract  lines  thru  keyboard  Instructions.  A 
subroutine  containing  the  equations  of  the  transfer  function  of  the  spectro- 
scope then  plots  the  effects  of  this  simulated  data  against  the  background  of 
the  actual  data.  When  the  user  is  satisfied  with  the  goodness  of  fit,  cards 
are  punched  out  containing  the  wavelength  and  Intensity  of  the  lines  as  de- 
rived from  this  particular  rocket  flight.  The  user  can  then  compare  this 
flight  with  other  flights,  with  laboratory  experiments,  or  with  the  results 
of  other  workers  in  the  field. 
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STUDY  OF  VARIATIONAL  PRINCIPLES  IN  APPROXIMATIONS 
OF  MOLECULAR  CHARACTERISTICS 


Initiator  : K.  Yoshino. 

Problem  No.  t 4718 

Project  No.  : 8627 


This  problem  la  Involved  with  the  study  of  aome  of  the  characteristics  of 
various  molecules  as  they  change  from  one  energy  level  to  another  as 
examined  from  a spectroscopic  discipline.  Since  a good  deal  of  perturbation 
la  involved  the  analysis  Involved  deala  with  the  basic  variational  principles 
to  best  achieve  tn>imir<l"**<ni1*  of  the  molecular  characteristics.  What  is 
ultimately  desired  is  the  evaluation  of  a particular  parameter,  H (see  below), 
which  minimizes  the  perturbation  In  the  energy  changes  so  that  calculation 
of  molecular  constants  which  represent  certain  natural  elements  at  various 
energy  levels  may  be  accurately  determined. 

A computer  program  was  produced  to  aid  in  determining  those  constants 
which  establish  at  what  energy  state  the  molecules  under  consideration  attain. 
The  results  of  the  approximation  analysis  were  then  compared  to  experimental 
spectroscopic  data. 

The  analysis  involved  is  basically  the  application  of  a linear  "least  squares” 
routine  which  employs  variational  principles,  i.  e. , a bracketing  type  search 
performed  on  a particular  parameter  which  minimizes  the  standard  deviation 
of  observed  data  to  calculated  data. 

For  each  value  of  this  variational  parameter  a "least  squares”  fit  is  done 
and  a standard  deviation  is  calculated  for  that  particular  value  of  the  parameter 
so  eventually  that  particular  value  of  the  variational  parameter  for  which 


- 29 


■■  ■ 


the  standard  deviation  is  minimized  is  then  taken  to  be  that  value  which  renders 
the  least  amount  of  perturbation. 

The  following  is  the  mathematical  sequence  involved  which  describes  the 
pertinent  spectroscopic  relationships  required  to  achieve  the  least  amount 
of  perturbation  as  various  energy  levels  are  achieved: 

1)  Tj(x)  = Aj  + BjX  - njX2 

2)  T2(x)  = Aj  + B2x  - D2x2 

3)  x = J(J  + 1)  J's  are  integers 

4)  WJ(x)  = (T^x)  + T2(x))/2  + (l/2)|/fi2  + 4H2  + when  6 > 0 

5)  W2(x)  = (Tx(x)  + T2(x))/2  ± (l/2)|/62  + 4H2  + when  6 < 0 

6)  6 = - T2 

7)  62  = (W*  - W*)2  - 4H2 

8)  T*(x)  = (W*  + W*)/2  ± 6/2 

9)  T*(x)  = (W*  + W*)/2  ± 6/2 

(a)  Find  6(J)'s  by  Eq.  (7)  for  OiHs  100. 

(b)  Find  T^(x)  and  T2(x)  by  Eq.  (8)  and  Eq.  (9).  Calculate  constants  A,  B, 
and  D in  Eq.  (1)  and  Eq.  (2)  using  values  of  T^(x)  and  T2(x)  obtained  in  (a) 

(c)  Calculate  the  standard  deviation  of  T^(x)  and  T2(x). 

(d)  Back  to  (a)  and  repeat  the  steps  to  find  that  H by  which  the  standard 
deviation  is  minimum. 

(e)  Calculate  T^x)  by  Eq.  (1),  Tg(x)  by  Eq.  (2),  and  then  calculate  W^x) 
by  Eq.  (4)  and  W2(x)  by  Eq.  (5)  using  the  value  of  H which  minimizes 
the  standard  deviation. 

Given  data:  W?(x),  Wi(x),  D, . and  D„. 


RADIO  ANTENNA  PROGRAM  CONVERSION 


Initiator  : Maj.  A.  Martinez 

Problem  No.  : 4719 

Project  No.  : 5635 


The  work  to  be  done  consisted  of  converting  a large  numerical  analysis 
program  for  radio  antennas  to  run  cm  the  CDC  6600. 

The  initiator  supplied  two  versions  of  the  program;  one  was  set  up  to 
run  on  an  IBM  360/65;  the  other  supposedly  could  be  run  on  the  CDC  6600, 
and  the  hope  was  to  reduce  its  core  requirements  and  increase  its  efficiency. 

The  CDC  version  turned  out  to  be  in  the  wrong  format  and  have  irretriev- 
able errors  in  the  character  set. 

The  IBM  version  was  treated  to  remove  syntactical  errors  and  remove 
incompatibilities  between  the  machines.  It  then  became  obvious  that  it 
would  be  easier  to  work  on  the  CDC  version,  line-by-line. 

Hundreds  of  corrections  were  made.  Despite  tremendous  core  memory 
requirements,  a few  simple  cases  were  run. 

In  the  final  analysis,  the  initiator  was  supplied  with  an  item-by-item  list 
of  why  the  original  tape  could  not  have  been  as  advertised,  what  coding  was 
known  bad,  and  what  was  suspect. 

In  spite  of  a large  amount  of  effort,  no  final  program  was  able  to  be 
submitted  to  the  SUYA  library. 
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MODEL  FOR  THE  REFRACTIVE  INDEX 
OF  STRATOSPHERIC  AEROSOLS 


Initiator 
Problem  No. 
Project  No. 


E.P.  Shettle. 

4731 

7621 


PURPOSE:  To  develop  an  analytic  expression  for  the  refractive  index 

of  the  stratospheric  aerosol.  This  permits  a more  compact  storage  of  the 
available  data  and  provides  an  accurate  extrapolation  to  wavelengths  for  which 
data  is  not  available.  The  refractive  indices  as  a function  of  wavelength  will 
later  be  used  to  determine  other  optical  properties  of  model  stratospheres. 
The  stratospheric  aerosol  is  predominately  a 75%  sulfuric  acid  solution. 


DESCRIPTION:  (1)  Do  a least  squares  fit  of  given  data  for  the  real 

and  Imaginary  parts  (n  and  k,  respectively)  of  75%  sulfuric  acid  solution 
versus  frequency  (v)  to  the  following  analytic  expressions: 


where  the  Aj,  Bj,  y.,  and  ia  (j  = 1, ...,  N)  are  parameters  whose  values  were 
to  be  determined  by  the  least  squares  optimization. 


For  approximate  values  of  the  parameters,  the  following  considerations 
were  used: 


the  t/j  fell  at  the  peaks  of  k(v) 

the  y^  are  the  width  of  die  peak  or  the  distance  between  the  corresponding 
turning  points  of  n(p) 


B 


j 


■ys 


This  leads  to  the  following  approximate  values  for  the  parameters: 
(j  — ••• ) 7)« 


L 

a 

1 

l 

450 

70 

38,000 

10, 500 

2 

580 

70 

60, 000 

25, 000 

3 

900 

50 

81,000 

22,000 

4 

1,050 

70 

170, 000 

50, 000 

5 

1,160 

200 

450, 000 

150, 000 

6 

1,  720 

21C 

225, 000 

80, 000 

7 

2,  880 

800 

900, 000 

350,000 

These  values  were  used  as  a starting  point  for  a least  squares  iteration. 


The  optimization  was  achieved  using  subroutine  ZXMARQ  which  minimizes 
the  sum  of  squares  of  functions  using  a modified  Marquardt  algorithm.  The 
program  was  first  used  on  the  real  and  imaginary  parts  separately  in  order 
to  generate  good  initial  guesses  for  the  coefficients.  These  were  then  used 
as  inputs  to  the  full  blown  model,  hi  all,  29  coefficients  were  used  for  the 
fit;  8 A.’s  and  7 each  of  B^,  y.,  and  u.. 

(2)  After  finding  the  correct  parameter  values,  calcomp  plots  were 
made  comparing  the  fitted  analytic  values  of  n and  k with  the  experimental 
data  versus  both  frequency  and  wavelength. 


33 


Initiator  : U.  Lamm  era. 


Problem  No.  : 4733 

Projeot  No.  : 8682 


The  problem  requires  evaluation  of  a navigation  system,  to  be  used  as  an 
aircraft  independent  landing  monitor  (ILM)  under  multipath  interference  con- 
ditions. For  simplicity,  a two-dimensional  configuration  with  a single  image 
source  is  assumed.  The  type  of  smoothing  to  be  used  an  the  phase  angle 
velocity  of  the  received  signals  will  be  studied.  Another  study  area  involves 
the  effects  of  this  smoothing  on  the  overall  system  accuracy. 


The  figure  illustrates  a case  where  a CW  signal  transmitted  by  M is  re- 
ceived by  two  antennas,  S and  S^,  which  are  spaced  at  a distance  s and  which 
are  in  linear  motion.  After  a time  interval  At  or  2At  the  antenna  system 
moves  from  the  right  position  to  the  center  and  left  positions,  respectively. 


We  call  k • i1  the  total  accumulated  phase  shift  measured  at  S,  when  S is  moving 
from  the  right  to  the  center  position  and  k • the  difference  in  total  accumulated 
phase  between  S and  along  the  same  path. 
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Correspondingly,  quantities  k • i g and  k • 6g  are  measured  when  moving  S 
and  Sj  from  the  center  to  the  left  position.  Here,  k is  2 ff/X  with  X the  trans- 
mitted wavelength.  If  X,  s,  At,  ip  ig,  6p  and  6g  are  known  the  other 
parameters  can  be  computed 


. + v 

lr212r^iT 


hl-«= 


WlV8!  + 62><‘l  * ‘2> 

S(62‘i  - S^)2 


v“f/1<i2-il>*°-5(1l+12)i 


Multipath  interference  on  this  system  is  simulated  by  an  image  source  M'  of 
amplitude  R.  With  signals 


Vs  = e*MS  + Re*5^  and 


jkMS-  jkM7S1 
Vc  = e + Re 


arriving  at  antennas  S and  the  phase  angular  velocities  of  the  signals  are 

= Re (VI  « d Im(Vl/dt  - Im(V)  • d Re(V)/dt 
Re2(V)  + Im2(V) 

where  Re(V)  and  Im(V)  represent  the  real  and  imaginary  part  of  the  Vg, 

respectively,  vL  . 

1 

Quantities  ip  ig,  and  6g  are  computed  through  numerical  integration 
of  ckp/dt.  Various  types  of  smoothing  of  dp/dt  under  a given  multipath  situa- 
tion are  evaluated  by  the  effect  that  this  smoothing  has  on  the  accuracy  of 
the  output  values  hj  - e,  jt,  and  v.  The  distance  l is  large  compared  with  s. 
In  this  case  rays  arriving  at  S and  may  be  considered  to  be  parallel. 


Convergence  problems  which  have  been  observed  are  believed  to  be  due 
to  the  basic  assumption  of  parallel  rays.  Here,  the  cases  of  no  multipath 


- 35  - 


and  perfect  reflectivity  (i.  e. , perfect  symmetry)  do  not  create  a problem. 
However,  when  partial  reflection  is  allowed  the  assumption  of  parallel  rays 
did  introduce  an  instability  in  the  altitude  estimate. 


WITHDRAWAL  WEIGHTING  SYNTHESIS 


Initiator  : 
Problem  No.  : 
Project  No.  : 


A.  Slobodnick 
4759 

634B-01-99 


We  are  concerned  here  with  the  generation  of  time  functions  for  withdrawal 
weighting.  That  is,  we  desire  the  overall  transfer  function  given  by  D(f), 
where  D(f)  is  realized  by  cascading  two  filters  Hw(f)  and  H^(f).  Thus, 


D(f)  = HA(f)  Hw(f)  . 


Now,  Hw(f)  is  fixed  and  specified.  Therefore,  we  must  realize  HA(f),  where 
HA(f)  is  given  by 


HA(f)“ir§) 


W' 


Now,  D(f)  is  to  be  an  ideal  bandpass  filter,  that  is, 


D(f)  = 1 within  the  given  band 
= 0 otherwise. 


Therefore,  we  are  attempting  to  realize 

H .(f)  = „ * within  the  bandwidth 
A'  Hw(f) 

= o otherwise  . 

We  approximate  HA(f)  by  a non-recursive  digital  filter  with  a frequency 
response  given  by 
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The  results  of  the  design  were  presented  to  the  investigator.  The  results 
were  considered  to  be  good  and  thus,  the  program  was  reformated  into  the 
form  desired  by  the  investigator.  The  final  program  was  given  to  the 
investigator  for  his  use. 


I 


MEASUREMENT  OF  AMBIENT  ELECTRIC  FIELD 

Initiator  : M.  Smiddy. 

Problem  No.  : 4764 

Project  No.  : 8617 


The  measurement  of  ambient  electric  fields  on  two  rockets  is  part  of  a 
large  AFCRL  program  involving  chemical  release  and  mass-spectrometer 
measurements. 

The  processing  involved  can  be  divided  into  three  phases  where  one  phase 
must  be  completed  before  commencement  of  the  next. 

I.  Digitizing  of:  1)  Rocket  trajectory  data. 

2)  Rocket  attitude  data. 

3)  Experiment  data  AV. 

All  input  data  will  be  on  magnetic  tape  and  computer  (CDC  6600)  output  tapes 
are  required. 

n.  Trajectory  and  attitude  data  to  be  processed  to  give  probe  position  and 
attitude  data  as  a function  of  time-from-launch. 

One  vehicle  (A10. 302-1)  has  three-axis  fluxgate  magnetometers  for  attitude. 
The  second  vehicle  (A10. 403-3)  has  magnetometers  and  a gyro  for  attitude 
acquisition  - this  flight  has  also  a mass-spectrometer  aboard  who  will  also 
require  attitude  processing. 

Ill  a.  Two  dipole  E-field  data  reduced  to  give  three  spatial  components 
Ey,  Ey,  and  Ez  also  components  along  and  perpendicular  to 
geomagnetic  field  B.  This  is  performed  in  a way  which  is  depen- 
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dent  on  the  vehicle  motion  but  ia  performed  by  subtracting  the 
Lorenze  field  produced  by  vehicle  motion  from  the  measured 
field  at  a time  t (above  SO  km  altitude). 


ft  • ft  = 

where 


-|Jy-<Vxfi).ft  , 

ft  - ambient  electric  field  vector 

ft  * dipole  position  vector 

|R|  = length  of  dipole 

V = voltage  differences  measured 

v * vehicle  velocity  vector 

ft  - geomagnetic  field  vector  . 


Other  diagnostic  Information  is  also  computed  such  as:  a function  which 
shows  the  magnitude  of  the  magnetic  field  connection  area  between  the  vehicle 
and  the  dipole  sensors;  the  velocity  veetor  location  in  relation  to  the  dipole 
axes. 


m b.  The  A10. 302-1  rocket  has  also  a spherical  electrostatic  analyzer 
for  measuring  electron  density  and  temperature.  The  equations 
«nH  mode  of  processing  this  data  will  be  made  available  as  re- 
quired (the  method  is  similar  to  data  processed  on  ISIS- 1 and 
INJUN-V  satellites). 


- 40 


CONTOUR  OF  EARTH  DEFORMATION 
DUE  TO  CHEMICAL  EXPLOSION 


Initiator  : M.  Settle. 

Problem  No.  : 4774 

Project  No.  : 7639 


Recent  research  on  deformation  in  granular  earth  materials  (e.  g. , soils 
ftivi  sands)  produced  by  chemical  explosions  is  being  conducted  within  the 
Terrestrial  Sciences  Laboratory.  Field  tests  conducted  nearby  at  Ft.  Devens 
have  been  designed  to  specifically  measure  in  situ  ground  strengths  both  be- 
fore and  after  the  detonation  of  a chemical  explosive.  Ground  strength 
measurements  are  recorded  as  continuous  'profiles'  on  a penetrometer  device 
which  probes  to  depths  up  to  20  inches  along  a single  profile.  The  purpose  of 
this  research  is  to  ultimately  delineate  the  extent  and  geometry  of  the  zone  of 
permanent  deformation  surrounding  the  crater  produced  by  such  explosions. 

Two  programs  were  written  to  do  the  above  described  job.  The  first  pro- 
gram generated  the  probe  curve  equations  of  ground  "strength*  values  for  a 
given  set  of  probes  and  combines  deep  input  probes.  It  consists  in  producing 
mathematical  equations  which  accurately  describe  ground  "strength"  measure- 
ments recorded  as  continuous  "profile"  on  a penetrometer  device  which  probes 
into  the  ground. 

The  second  program  uses  the  output  of  the  first  to  generate  a contour  plot 
of  ground  "strength”  values  within  given  horizontal  and  vertical  distances 
from  the  center  of  explosion.  The  contour  plot  is  drawn  from  a set  of  given 
values  produced  from  the  input  probe  equations.. 
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LAMB  WAVE  VELOCITY  DISPERSION 


Initiator  : T.  L.  Szabo 

Problem  No.  s 4787 

Project  No.  : 5635 


i 


Lamb  waves  are  accoustic  waves  that  propagate  in  plates  at  a velocity 
determined  by  its  frequency  and  the  thickness,  density,  and  elastic  constants 
of  the  plate.  Calculations  of  the  velocity  of  the  first-order  symmetric  and 
and -symmetric  Lamb  modes  are  required  for  nondestructive  evaluation  of 
stressed  plates. 

For  a given  type  of  elastic  solid  the  longitudinal  (c^)  and  bulk  (c^)  wave 
velocities  are  known.  The  velocities  of  the  first-order  symmetric  (Bq)  and 
anti-symmetric  (aQ)  modes  are  to  be  calculated  and  plotted  as  a function  of 
the  running  parameter  2fh,  which  is  the  frequency  (f),  thickness  (2h)  product 
of  the  plate  in  which  these  Lamb  waves  propagate. 

The  following  equations  govern  the  velocity: 

O* 

If  u =-jr  and  V = 

CZ 

the  first-order  symmetric  mode  equation  is 

tan  (1  - u)l/2  ■£-  (2fh)  j/2  1/2 

S-  4u(l  - u) ' (v  - *)  ' 

tan  (v  - u)1/2  (2fh)  (2u  - l)2 

and  the  first-order  antisymmetric  mode  equation  is 

tan  (1  - u)l/2  7?-  (2fh)  2 

W (2u  - 1) 

tan  (v  - u)1/2  (2fh)  4u(l  - u)1/2  (v  - u)1//z  * 
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These  equations  are  solved  for  c using  the  root  finding  routine  ZREAL2. 


The  program  has  been  run  for  a variety  of  materials  and  parameters  as 
requested  by  the  initiator. 


ACOUSTIC  WAVE  BEAMSHAPING 


Initiator 
Problem  No. 
Project  No. 


T.L.  Szabo. 

4788 

5635 


This  problem  deals  with  acoustic  wave  beamshaping  for  non-destructive 
testing.  (Ultrasonic  non-destructive  evaluation  of  defects  is  complicated  by 
changes  in  the  beam  pattern  caused  by  diffraction. ) 

18  versions  of  programs  were  written.  10  of  these  are  in  the  process  of 
being  submitted  to  the  SUYA  library. 

These  computer  programs  deal  with  aspects  of  simulating  various  wave- 
shapes as  they  propagate  along  a material.  They  progress  to  one  which  can 
analyze  any  arbitrary  shape  at  the  launching  transducer  and  culminate  in 
one  which  reconstruct  what  the  signal  must  have  looked  like  at  the  launching 
transducer,  from  the  way  it  looked  when  picked  up  further  along  the  material. 

The  goals  originally  envisioned  in  the  statement  of  work  have  more  than 
been  accomplished.  Any  of  the  various  encouraging  possibilities  opened  up 
by  these  programs  will  have  to  be  re-initiated  as  a separate  problem  in  the 
future. 
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A CALCULATION  OF  FOURIER  HARMONICS 


Initiator  : A.  Cole 

Problem  No.  : 4793 

Project  No.  : 8624 


This  problem  is  involved  in  calculating  AQ,  A^,  Ag,  <f>^,  and  <P2  in  the 
following  equation: 

y(x)  = AQ  + Aj  sin  (x  + <p±)  + Ag  sin  (2x  + <P2) 


where 


x = iZ,  Z = 360/period 
i = 0,3,6,9,12,15,18,21 
Aq  = mean 


A1 

*2 


= amplitude  of  the  first  harmonic 
= amplitude  of  the  second  harmonic  . 


A table  of  values  of  y's  are  to  be  supplied  by  the  initiator. 


A Fourier  analysis  computer  program  was  developed  to  successfully  pro- 
duce the  desired  harmonics.  The  results  of  this  effort  were  used  in  the 
overall  patterns  of  a meteorological  study. 


INTEGRAL  COMPUTATIONS  USING  GAUSS-MEHLER  QUADRATURE 


Initiator 
Problem  No. 
Project  No. 


D.  Freeman 

4808 

6687 


PROBLEM  DESCRIPTION  Experimental  values  are  given  for  the  discrete 
vibrational  energies  G(v)  with  v = 0, 1, ... , v and  for  the  rotational  con- 

III  HX 

stants  B(v)  with  v = 0, 1, ...,vmax.  We  wish  to  compute  the  integrals: 


f = 4. 1057886M"1/2 

V 

f 

[G(v)  - G(v')f 1/2  dv' 

and 

min 

gv  = 0.24355857M1/2 

f 

B(v/)[G(v)  - G(v')fl/2 

✓v 

vmin 

using  Gauss-Mehler  quadrature,  as  described  by  Tellinguisen.1  This  tech- 
nique is  fast  and  avoids  the  singularity  at  the  upper  limits  of  integration. 
The  reduced  mass  M of  the  oscillator  and  the  value  of  vmnv  are  given.  The 
lower  limit  of  integration  is 


v = _ A _ 6 
min  2 

with 

6 = (2x)-1[(w2  + 4xY)l/2-  w) 

The  values  of  G(v)  and  B(v)  at  v = vmin  are,  if  6 ^ 0, 


G<W  = ’Y 
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(7) 


= B + g6  + y62  . 

In  Eqs.  (4) — (7),  the  vibrational  parameters  w and  x,  and  the  rotational  con- 
stants B,  a,  and  y are  given.  If  6 = 0,  set 

G<W=Y  = 0>  B<W * B • 

The  program  provides  the  option  of  constructing  G(v')  and  B<v/)  in  Eqs. 
(1)  and  (2)  in  either  one  or  both  of  the  following  ways: 

(a)  G(v')  and  B(v')  are  constructed  as  piecewise  analytic  functions  of  v'  by 
performing  cubic  Hermite  interpolation  between  consecutive  experi- 
mental values,  and  the  required  derivatives  at  each  such  point  are 

o 

formed  artificially  by  a finite  difference  approximation. 

(b)  G(v')  and  B(v')  are  represented  by  the  power  series 

6 , Nl-1 

B<v'>=£Crl(','+§) 
i=l 

in  which  the  coefficients  are  supplied  as  input. 

It  is  known  that  G(v)  increases  monotonically  with  v and  usually  has  a 
negative  second  derivative  and  that  B(v)  decreases  monotonically  with  v and 
usually  has  a negative  first  derivative.  The  coefficients  in  Eqs.  (8)  and  (9) 
decrease  rapidly  in  magnitude  as  i increases. 

OUTPUT  A title  such  as  "RKR  curve  for  the  ground  state  of  CO" 

Reduced  mass  M = • • • • 

Vibrational  parameters: 


(8) 

0) 
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a listing  of  the  coefficients  in  Eq.  (8)  if  used 
Rational  constants: 

a listing  of  the  coefficients  tn  Eq.  (9)  if  used 
A tabulation  of  the  following  quantities: 

v,  B(v),  G (v)  + Y,  2fy,  20gy,  Ry  (outer),  Ry(inner)  for  the 
values  v = ne  with  n = 0, 1, 2, ... , vmnT|/  A step  size 

c is  given,  normally  in  the  range  0.25  - 5.0.  The  quantities 
R (outer)  and  R (inner)  are  defined  by 


+ f and  Ry  (inner) 


Rfouter) 


Also,  given  as  output  are  the  quantities:  6,  Y,  R(vmin) 

quantity  R(v_.  ) is  defined  as 


4. 1057886M 


SAMPLE  PROBLEM  Data  and  results  for  the  ground  state  of  CO  are  con- 
tained in  the  references. 


REFERENCES  (attached) 


1.  J.  Tellinguisen,  Computer  Phys.  Commun.  6,  221  (1974); 

J.  Tellinguisen,  J.  Mol.  Spectrosc.  44,  194  (1972). 

2.  P.  Tsipouras  and  R.  V.  Cormier,  AFCRL-TR-73-0400  (1973) 

3.  A.W.  Mantz,  et  al..  J.  Mol.  Spectrosc.  39,  180  (1971). 


COMPUTER  MODIFICATION  USED  IN  THE  CALCULATION  OF 
RADIAL  AND  TANGENTIAL  DISPLACEMENTS 
OF  THE  ELASTODYNAMIC  FIELD 


Initiator 
Problem  No. 
Project  No. 


P.  Michaels 

4813 

7639 


This  problem  deals  with  the  modification  of  a computer  program,  ELAS, 
which  calculates  the  radial  and  tangential  displacements  of  the  Elastodynamic 
field  radiated  by  an  impulsive  force  acting  normal  to  a line  in  the  surface  of 
a semi-infinite  medium. 

This  problem  is  a continuation  of  a problem  that  is  described  in  detail  in 
AFCRL-68-0481.1  The  total  problem  is  concerned  with  the  detection  and 
analysis  of  underground  seismic  waves. 

The  determination  of  stress  components  of  a cylindrical  S-wave  and  P- 
wave  propagated  in  a pre-stressed  medium  has  been  determined.  The 
investigation  and  analysis  has  been  completed  in  the  determination  of  stress 
components  which  involve  the  following  simulated  stresses: 

• a stress  field  around  a hole  in  an  otherwise  unaxial  compression  field  (1), 

• a pure  shear  field  around  a Held  (2), 

• a pure  tension  field  around  a hole  (3), 

• and  any  additive  combination  of  1 thru  3. 


Richard  W.  Doherty,  Joseph  F.  Russell,  and  Rose  M.  Ring,  Scientific 
Research  in  the  Form  of  Numerical  Analysis  and  Data  Analysis.  AFCRL 
68-0481,  Bedford,  Mass.,  p.  30. 
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Another  force  was  introduced  into  the  pre-stressed  medium.  This  new 
force  was  a dipole  radiation.  Hence,  a program  was  written  to  determine 
stress  components  for  dipole  radiation.  These  components  were  determined 
from  the  following  equations: 

or  . o(.  - *)  - ApCpi  »(.-$) 

ad  = bar 

TrS,Adc.22sJ»Gl(t-£)  , 

where:  If  0 £ [t  - (r/Cp)]  £ then, 
otherwise,  F(t  - (r/Cp)l  = 0. 


If  [t  - (r/Cp))  a 0,  then 


otherwise,  G[t  - (r/Cp)]  = 0. 

If  [t  - (r/Cs)]  a 0,  then 


otherwise,  G[t  - (r/Cs)]  = 0. 

I. 


Tho  Mnim  rtf  the  narnmpforfl  la  an  fnllmvat 


. 


0®  s 0 * 90®  A0  » 10® 

1.0  *r<  2.0  A(r)  = .1 

2.0  <r«90  A(r)  = 1.0 

Ad  - .5,  Ap  = 1.0,  T = 1.0,  and  tQ  = .4,  6 = .4 

The  required  modification*  were  made  to  the  computer  program,  ELAS, 
mentioned  above  nTtf1  was  subsequently  submitted  to  the  Computer  User's 
Library  at  AFGL,  Hanscom  Field. 


TAPE  EVALUATION 


Initiator  L.  Logan 

Problem  No.  t 4815  , 

Project  No.  : 7630 



The  problem  consisted  of  evaluating  a tape  produced  by  a specialized 
balloon  recording  system,  acceptance  of  the  hardware  by  the  Air 

Force. 

Several  programs  were  written,  and  various  tape  dumps  were  made 
at  Mitre.  Unfortunately,  a total  of  6 tapes  were  finally  supplied,  none 
of  them  being  in  the  advertised  hardware  format. 

Thus,  the  original  goal  of  evaluating  the  consistency  of  the  data  be- 
came subordinate  to  getting  any  proper  data  at  all. 


The  latest  program  was  HsnrioH  over  to  the  initiator,  for  him  to  incor- 
porate in  any  future  testing  efforts.  No  program  was  submitted  to  the 
SUYA  library  since  there  never  was  any  proper  data  supplied. 


POTENTIAL  BETWEEN  TWO  RARE  GAS  ATOMS 


i 


Initiator  D.  Freeman 

Problem  No.  : 4817 

Project  No.  : 6687 


PROBLEM  DESCRIPTION  The  distribution  of  the  eigenvalues  E(v)  of  a 
potential  V(R)  is  given  fay  (Ref.  1) 


dv 

dE(v) 


4 


^(v) 


(v) 


(E(v)  - V(R)fl/2  dR  , 


(1) 


in  which  A is  a known  constant  and  E(v)  is  the  eigenvalue  of  the  discrete 
level  v which  has  inner  and  outer  turning  points  R^(v)  and  Rgtv),  respec- 
tively. A simple  approximation  to  be  long-range  potential  between  two 
rare  gas  atoms  is 

V(R)  = D - CgR-6  , (2) 

in  which  D is  the  unknown  dissociation  limit  and  Cg  is  a known  constant 
related  to  the  outer  turning  point  by  the  expression 

E(v)  = D - CgR“6  . (3) 

Reference  (1)  shows  that  the  substitution  of  Eqs.  (2)  and  (3)  into  Eq.  (1) 
leads  to 

afe> * ac61/2r5  ( 2/  1 »'V  - 1>'1/2  * , w 
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• in  terms 


in  which  y = RgR-1  and  the  integral  Is  known  in  the  limit  Rg/Rj^  -* 
of  gamma  functions.  Equations  (3)  and  (4)  then  yield 


(5) 


in  which  Kg  is  a collection  of  known  constants.  If  v = vQ  is  defined  as  the 
vibrational  index  (generally  nonintegral)  at  the  dissociation  limit,  E(Vp)  = D, 
the  integrated  form  of  Eq.  (5)  is 


D - E(v)  = 


(6) 


The  problem  is  to  find  analogous,  preferably  in  analytical  form,  of  Eqs. 
(4),  (5),  and  (6)  when,  instead  of  using  the  pair  of  Eqs.  (2)  and  (3),  we  use 
<me  of  die  more  accurate  pairs 


V(R)  = D - C6R"6(1  + 0R~2) 

(7a) 

E(v)  =D-C6I^6(1  + 0RJ2) 

(7b) 

V(R)  = D - CgR-6  exp  (/5R“2) 

(8a) 

E(R)  * D - CgR"6  exp  (0Rj2) 

(*) 

V(R)  * D - C6R‘6(1  - fiR'2)’1 

<8a) 

E(R)  =D-C6R^6(1  -0RJ2)"1  , 

(8b) 

in  which  0R~2  and  0R^2  are  small  ( < 0. 2)  compared  to  imity.  We  require 
the  problem  to  be  solved  only  for  any  one  of  Eqs.  (7),  (8),  or  (9),  not 
necessarily  for  all  three. 


POLYNOMIAL  FIT  A polynomial  approximation  of  the  form 

V(R)  * D - oR"V  + PR~2  ♦ yR"4) 

E(v)  = D - cR^6(l  +0RJ2  ♦yRj4)  , 
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in  adopted  for  integral  (1)  by  which  it  becomes  a function  of  P and  y.  A bi- 
variate polynomial  is  fit  to  the  computed  values  of  dv/dE(v). 

REFERENCES 


1.  R.J.  Le  Roy  and  P.B,  Bernsten,  J.  Chem.  Phy.  52,  3869  (1970). 

2.  S.M.  Kirschner  and  J.  K.  Watson,  J.  Mol.  Spectroscopy,  47,  2346  (1973). 


NUMERICAL  MAPPING  OF  IONOSPHERIC  PARAMETERS 
ON  A GLOBAL  BASIS 


Initiator  : J.  Kotelly 

Problem  No.  : 4825 

Project  No.  : 7663 


The  problem  is  to  investigate  methods  of  numerical  mapping  ionospheric 
parameters  on  a global  basis;  specifically,  this  involves  developing  the 
capability  to  plot  values  of  FOF  2 as  a function  of  geomagnetic  latitude  for  a 
constant  local  mean  time. 

This  problem  was  acquired  only  a few  weeks  before  the  end  of  the  period 
covered  by  this  report,  and  will  be  carried  into  a future  one. 

Several  computer  programs  were  written  for  the  set  of  data  presented, 
but  they  only  proved  that  the  data  had  many  transcription  errors. 

Several  computer  programs  were  written  for  a second  set  of  data,  but 
progress  must  be  halted  until  the  cards  designating  the  observing  stations 
arrive  at  AFGL. 


TAPE  REFORMATTING 


Initiator  : E.  Robinson 

Problem  No. : 4826 

Project  No.  : 6690 


This  problem  Involves  reformatting  S AMT  EC  TRAE  radar  tracking  data 
from  IBM  7094  binary  tape  into  the  standard  SUYA  CDC  6600  radar  tape 
format. 

A program  was  written  to  perform  the  above,  and  to  produce  time  plots 
of  azimuth,  altitude,  and  range  as  well. 

Continued  minor  changes  have  been  made,  mostly  to  augment  the  control 
the  user  has  over  what  gets  plotted,  since  in  practice  it  seems  the  radar 
track  spurious  objects  before  and  after  the  desired  flights. 


INVERSION  OF  LIMB  AIRGLOW  MEASUREMENTS 


Initiator 
Problem  No. 
Project  No. 


R.  A.  Van  Tassel. 
4828 
6688 


This  analysis  provides  support  for  an  optical  mission  on  the  Defense 
Meteorological  Satellite  System  (DMSS).  The  objectives  are: 

1.  Computation  of  theoretical  optical  volume  emission  rates  as  a Amotion 
of  altitude; 

2.  Computation  of  the  simulated  output  of  the  DMCS  sensor  system;  and 

3.  Development  of  an  inversion  algorithm  to  infer  ionospheric  parameters 
from  DMSS  data. 

PART  I.  Theoretical  Optical  Volume  Emission  Rates 

The  computation  of  theoretical  optical  volume  emission  rates  as  a function 
of  altitude  will  be  computed  for  two  different  mechanisms:  1)  radiative  re- 
combination, and  2)  electron  impact  excitation. 

For  radiative  recombination,  we  will  use  the  formula 

F(z)  = orn(z)2  , 

where  F(z)  is  the  volume  emission  rate,  a is  the  rate  coefficient,  and  n(z)  is 

3 

the  number  of  electrons/cm  , given  by  the  relation 

In  this  relation,  S is  a shape  factor  and  He  is  the  electron  scale  height. 
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Initial  computations  may  use  values  as  follows: 

a = 8 x 10-13  cm3  sec”1  , 

5 -3 

n = 10  electrons  cm  , 

max 

zmax  = 250  km> 

SH  = 50  km  . 

6 

The  required  range  of  F(z)  is  from  200  to  500  km  in  increments  of  1 km. 

-3 

The  units  of  F(z)  are  photons  cm  . 

In  order  to  compute  theoretical  emission  rates  resulting  from  electron 
impact,  we  use  the  formula, 

F(z)  = nQ(z)  J <r(f)  ne(z,e)  dc  , 

-3 

where  n.(z)  is  the  number  density  of  atomic  oxygen  in  cm  , o(c)  is  the  cross 

^ 5 0 

section  for  electron-impact  excitation  of  the  S state  in  atomic  oxygen  as  a 

function  of  electron  energy,  and  ne(z,  c)  is  the  electron  energy  distribution 

function.  The  range  and  increments  are  again  200  to  500  km  in  1 km 

increments. 

The  function  nQ(z)  is  available  from  a model  atmosphere  which  is  already 
stored  in  the  computer.  An  experimental  determination  of  a(c)  is  shown  in 
Figure  1.  The  electron  energy  distribution  is  being  calculated  by  another 
program  in  this  work  unit.  Paul  Tsipouras  (SUYA)  is  familiar  with  this 
program. 

PART  II.  Simulated  Satellite  Output 


The  computation  of  the  simulated  output  of  the  DMSS  sensor  system  involves 
computing  the  column  emission  rate  by  integrating  along  lines-of-sight  or  ray 


paths.  The  effects  of  absorption  along  the  ray  path  is  included.  The  ex 
pression  for  the  column  emission  rate  4trt(h)  is  given  by 


4irl(h)  = 4irt+(h)  + 4irl_(h)  , 


where 


and 


4»I±(h)  =jf 

/■dx\  (Re  * z> 
\dz/  A(h,  z) 


• -WM)/dx 
F(z)  e A — 


(if)  <**  - 


with  A(H,  z)  given  by 


A(h,  z)  = ^ - h^  + 2Rg(z  - h)  , 

where  h identifies  the  ray  path  by  its  tangent  altitude  and  R is  the  radius  of 
the  earth. 

Initially  we  will  assume  an  isothermal  atmosphere  above  200  km,  so  that 
we  may  express  N^(h,  z),  the  number  of  absorbers  along  the  path  from  the 
emitting  point  to  the  sensor  located  mathematically  at  infinity,  by 

N+(h,  z)  = n0^-A(h,  z)  e“z/H  + (1)  J A(h,  $)  e_4/H  d^ 

N_(h,  z)  = n0^-A(h,  z)  e“z/H  + (-0  jf  A(h,  $)  e_€/H  d£ 

+ (I)j\(h,Oe^H<i{J  , 

where  H is  the  scale  height  for  molecular  oxygen. 
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These  results  should  be  compared  with  results  using  the  expressions 


f"  (R  + €) 

NA*>=J 

and  z 

f*  (Re  + £)  fh  (Re  + 0 

n(C)‘A(h^"d?+Jz  n(€)_A(h^rde  ’ 

where  n(£)  is  the  Og  number  density  from  a model  atmosphere  which  is  stored 
in  the  computer. 

PART  HI.  The  Development  of  an  Inversion  Algorithm 

The  first  two  parts  of  this  problem  were  concerned  with  calculating  theo- 
retical values  for  F(z)  and  then  4ffl(h)  from  a relationship  of  the  form 

4ffl(h)  = j F(z)  k(h,  z)  dz  . 

The  third  part  of  this  problem  is  the  inverse  of  the  second  part,  i.  e. , to 
infer  F(z)  given  several  observations  for  4frt(h).  The  general  approach  to  this 
kind  of  retrieval  problem  has  been  to  assume  an  initial  form  for  F(z)  and  to 
compare  calculated  values  of  4ffl(h)  with  observed  values.  Iteration  is  con- 
tinued until  the  desired  degree  of  accuracy  between  calculated  and  observed 
values  of  4irl(h)  is  obtained. 
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COLLIMATOR  PROGRAM  CONVERSION 


Initiator  : R.  Fredrickson. 

Problem  No.  : 4852 

Project  No.  : 5621 


This  job  was  essentially  a conversion  job.  A Fortran  program  written 
previously  to  run  on  another  computer  was  converted  to  run  on  the  CDC  6600 
computer  here.  This  program  gives  the  results  of  a Monte  Carlo  calculation 
of  the  energy  spectrum  and  angular  distribution  of  a gamma  radiation  emitted 
by  a collimated  Cobalt-60  source.  The  program  was  developed  flexible  enough 
to  give  results  for  a wide  variety  of  geometric  situations  of  source  container 
and  collimator. 

In  making  this  conversion  a Fortran  utility  routine  had  to  be  used  to  re- 
place an  existing  assembly  language  subroutine  which  was  untransferable. 

The  program  was  run  with  several  test  cases  varying  certain  input  param- 
eters to  obtain  several  sets  of  output.  After  close  examination  of  this  output 
the  initiator  was  assured  of  the  correctness  of  the  program. 


SIGNAL  ANALYSIS  OF  ROCKET  DATA 


Initiator  T.  Conley 

Problem  No.  : 4654 

Project  No.  : 7660 


This  problem  is  to  develop  a computerized  system,  including  appropriate 
computer  programs  and  operator  instructions,  to  digitally  filter  (smooth) 
unwanted  noise  from  rocket  measured  radiometer  and  photometer  data  (in 
digital  form),  correct  data  for  rocket  aspect,  and  digitally  differentiate 
smoothed  data  to  obtain  volume  emission  rates.  Current  ICECAP  rocket 
radiometer  measured  data  is  contaminated  by  noise  induced  by  the  telemetry 
recording  and  digitization  processes  and  which  leads  to  erroneous  volume 
emission  rates  when  differentiated. 

Considerable  data  exists  at  OPR  (approximately  10  rocket  flights)  which 
need  analysis  by  this  technique. 

Sample  data  was  copied  from  tapes  and  programs  were  developed  to  per- 
form initial  processing  on  the  data.  Specifically,  a noise  filter  and  differ- 
entiation filter  programs  were  written  and  used  on  the  data.  These  test 
results  were  in  excellent  agreement  with  known  results.  The  net  noise 
suppression  obtained  [noise  filter  minus  differentiation]  was  greater  than 
60  db. 

Thus,  development  of  a complete  system  is  underway  as  outlined  below. 

Work  on  this  problem  will  be  continued  under  another  contract. 
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Processing  Procedure  for  Emission  Rate  Determination 


Establish  data  file. 

Edit  data  (replace  calibration  and  noise  "spike”  paints  in  data) 
using  optimal  fitting. 

Convert  to  Brightness  Unit  (Ba). 

Plot  resulting  R a data  vs.  time. 

Produce  and  plot  FFT  of  data. 

Van  Rhijn  Aspect  Correction  (RQ  = Ra/i/a ). 

Plot  resulting  Ra  data  vs.  time. 

Produce  and  plot  FFT  of  data. 

Design  a noise  filter  and  process  data  through  the  filter 
(possibly  separate  filters  for  precession  and  wideband  noise). 
Plot  filtered  data  (vs.  time). 

Differentiate  output  of  noise  filter  (uBing  a differentiation  filter). 
Plot  resulting  output  (vs.  time). 

Convert  to  desired  variables  (emission  rate  vs.  altitude)  and 
produce  plot. 


i 


WAVELENGTH  PEAK  STUDY 


* Initiator  : T.  Conley 

Problem  No.  s 4855 
c Project  No.  : 7670 


This  problem  is  to  develop  a computerized  system  to  smooth  noise  con- 
taminated spectral  features  in  rocket  measured  CVF  data  (ICECAP  data)  in 
order  to  accurately  specify  wavelength  of  peak  emission  features  and  deter- 
mine line  widths.  This  technique  is  required  to  analyze  existing  data 
(approximately  5 rocket  flights). 

A noise  filter  is  being  developed  fox-  this  data.  However,  no  data  proc- 
essing has  been  begun  to  date.  Also,  techniques  for  "optimal*  estimation 
of  peak  locations  are  being  studied.  In  particular,  the  possible  use  of  full 
band  differentiations  (with  variance  estimates)  is  being  studied. 

Work  on  this  problem  will  be  continued  under  another  contract. 


STATISTICAL  ANALYSIS  OF  GAS  MEASUREMENT  DATA 


Initiator  : 
Problem  No. : 
Project  No.  s 


C.C.  Gallagher 

4861 

6687 


The  problem  is  to  develop  statistical  techniques  to  identify  gas  measure- 
ments from  column  and  spectrometer  tests.  Using  calibration  and  test  runs 
procedures  for  specifying  an  identification  procedure  is  to  be  developed  in- 
cluding probability  measures  of  the  uncertainties. 

A program  been  developed  to  calculate  the  "statistics1'  of  the  calibra- 
tion and  test  runs  provide  student  t tests  for  most  probable  compounds. 
Sample  data  will  be  used  to  evaluate  the  approach. 

Work  on  this  problem  will  be  continued  under  another  contract. 
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